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The equation of state and the shock Hugoniot of deuterium are calculated using a first-principles 
approach, for the conditions of the recent shock experiments. We use density-functional theory 
within a a classical mapping of the quantum fluids [Phys. Rev. Letters, 84, 959 (2000)]. The 
calculated Hugoniot is close to the Path-Integral Monte Carlo(PIMC) result. We also consider 
the quasi- equilibrium two-temperature case where the Deuterons are hotter than the electrons; the 
resulting quasi-equilibrium Hugoniot mimics the laser-shock data. The increased compressibility 
arises from hot D + — e pairs occurring close to the zero of the electron chemical potential. 
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Hydrogen and its isotopes have been extensively stud- 
ied, yet laset-shock experiments of l)a Silva et al.cl, 
Collins et alE, and Mostovych et al.B produced unex- 
pected disagreement with the equation of state (EOS) 
of the SESAME databaseo. The disagreement with 
SESAME occurs for temperatures T with ~ 0.8 eV 
< T < ~ 10 eV, and for densities 1.8 < r s <~ 2, where 
the electron-density parameter (in atomic units) is, r s — 
(3/4-Tm) 1 / 3 . Here n is the electron-number density (au.). 
The coupling constant T= (potential energy) /(kinetic en- 
ergy) ranges from 1 to ~ 30. 

The electrons change from a degenerate liquid to a clas- 
sical system in the anomaly-regime ( AR) , while the e — 
D + interaction is at the threshold of bound-state forma- 
tion. Thus the AR poses a difficult, strongly-correlated 
many-body problem of wide interest - from astrophysics 
and fusion to materials science. Hence a flurry of activity 
has focused on the deuterium EOS and its Hugoniot.El 
These involve intuitive approaches (called "chemical 
models"), til assuming the existence of molecules etc., and 
first-principles approaches, e.g., density- functional the- 
ory (DFT), or quantum Monte Carlo (QMC). Stan- 
dard DFT treats the ions as an external potential, opti- 
mized via Car-Parrincllo type techniques, while the local- 
density- approximation (LDA), or simpinj; tight-binding 
approaches are used for the electrons BeI Path-integral 
Monte Carlo (PIMC), a finite-T QMC approach, has 
been used for this problem!] 

PIMC and DFT calculations yield shock Hugoniots 
close to SESAME, showing no strong softening. Recent 
experiments from Sandia are also close to the SESAME 
Hugoniot .113 However a high degree of agreement exists 
in the laser-shock results.tll 

Our objective is to present new first-principles results 
for the EOS and Hugoniot of deuterium, using methods 
quite independent of previous methods. Our method is 
computationally simple and uses one-dimensional inte- 
gral equations based on DFT. We calculate all the pair- 
distribution functions (PDFs), gij, with i =1,2,3; (or ej, 
ej, and D + nuclei) of a three component system. Being 
based on the PDFs, it is manifestly non-local and has 



no self-interaction errors, ft easily adapts to the quasi- 
stationary two-temperature case with To 7^ T e . Consider 
the fluid with a deuterium nucleus at the origin, and let 
the one-body densities of the electrons and D + be n(r) 
and p{f). Then n(r) is really nr> e { r ) = ngr> e (r). The 
free energy F is a functional of the form F[n(r), p(r)]. 
Thus, taking functional derivatives, we have two coupled 
Kohn-Sham-Mermin equations: 

5F[n(r),p(r)]/5n(r) =0 (1) 
5F[n(r),p(r)]/5p(r)=0 (2) 

As shown in ref. [ |l2|;[l3) Eq.(l) leads to a (quan- 
tum) Kohn-Sham (K-S) equation for the electrons, while 
Eq.(2), a classical K=S equation, becomes the Hyper- 
Netted-Chain (HNC)Ell equation for a specific choice of 
the correlation potential (there being no exchange po- 
tential in the classical system). Thus the K-S eigen- 
functions, as well as the n(r), p(r), in thepHydrogen 
fluid were calculated by solving Eqs. (1-2). t 2 ] In Car- 
Parrinello approaches only the electrons are treated in 
DFT, while the N ions at sites Ri are explicitly treated. 
Here both electrons and ions are treated using their dis- 
tributions. Our hydrogen calculations were later con- 
firmed by lengthy QMC.E3 However, while Eqs. (1-2) pro- 
vide the gie( r ), 9ii(r), i=ion, and the K-S states in the 
fluid, the g e e(r) was available only in LDA. 

Recently, we showed how the electrons at the physi- 
cal temperature T e could be replaced by an equivalent 
classical system at T ee — (T 2 + T^) 1 / 2 , such that the 
quantum effects are correctly incorporated. A simple ex- 
pression for the electron "quantum temperature" T eq as 
a function of r s was given in ref. [l^. Application of the 
method, denoted the "Classical mapping of quantum sys- 
tems using the Hyper-Netted-Chain equation (CHNC)" , 
to 3-D and 2-D uniform electron liquids at T — and 
finite-T, showed excellent agreement of the g^, ener- 
gies, etc., with .QMC results, for even very strongly cou- 
pled situations JI3TL3 Deuterons in a uniform neutralizing 
background are mathematically identical to the electron 
system, except for changes of scale, (e.g, for deuterons, 
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r S D = r s (MD/m e ) etc., Mr> is the deuteron mass and 
m e — 1). Hence the -D + -quantum temperature Tjj q is 
also available, and is negligible in the regime of interest; 
thus D + are treated as classical particles. 

The mean densities p and n are equal since the nuclear 
charge Z = 1. Consider a fluid of total density n to t, with 
three species. Let Xi = iii/ntot, ntot — p + n. The physi- 
cal temperature is T, while the inverse temperature of the 
electrons is l/0 ee , with l/j3 ee = V(T 2 + T 2 ). For D+ , the 
classical mapping Tdd — 1/Pdd is taken to be T. The 
inverse temperature /3y between electrons and deuterons 
is actually not needed since the product, (3 e D<f>eD{i') is 
completely determined by the quantum mechanics of the 
problem (see below). However, for the present applica- 
tion it is the mean kinetic energy of the D + — e pair, 
i.e., T e o = (Tee + Tfjrj)/2, consistent with the classical 
picture used here. 

The classical equations for the PDFs and the Ornstein- 
Zcrnikc(OZ) relations are: 

9ij(r) = expl-Pijfajir) + hij(r) - (r) + By(r)] (3) 

hij(r) = cy(r) + Z s n s J dr'h i<s (\r - r'|)c SJ (r') (4) 

Here (j>ij{r) is the pair potential between the species i,j. 
For e-e (or D + — D + ) this is just the Coulomb poten- 
tial V cou (r). If the electron spins are parallel, the Pauli 
principle prevents occupation of the same spatial orbital. 
As beforejla we introduce a "Pauli exclusion potential" , 
V(r). Thus 4>ij(r) becomes V(r)Sij + V cou (r), when i, j 
denote electrons. The function h(r) = g(r) — 1 is related 
to the structure factor S(k) by a Fourier transform. The 
c(r) is the "direct correlation function (DCF)" of the OZ 
equations. The -By(r) is the "bridge" term due to cer- 
tain cluster interactions. If this is neglected, Eqs. [|- ^ 
form a closed set defining the HNC approximation. (In 
effect, the classical K-S equation becomes the HNC equa- 
tion if the correlation potential is evaluated as a sum of 
hyper-netted-chain diagrams.) The HNC is sufficient for 
the uniform 3DEG_fjor a range of r s , up to r s — 50, as 
studied previouslyO We neglect the bridge corrections 
in this study of deuterium. 

The V(r) is defined as in rcf. from the zeroth-order 
PDFs of the parallel-spin electrons. Thus: 



PV(r)=h 11 (r)-c 11 (r)-H9° 11 (r)} 



(5) 



where, e.g., c? 1 (r) is the spin-|| DCF of the O-Z equation. 
Only the product (3V(r) is needed. 

The Coulomb potential V cou (r) for two point-charge 
electrons is 1/r. However, an electron at the tempera- 
ture T is localized to within a thermal de Broglie length 
(dBL). Thus,pi)r the 3DEG we used a "diffraction cor- 
rected" formrj, 

V% u (r) = (l/r)[l-e- rk °°] (6) 

where k ee — {2itm*T ee ) 1 / 2 as shown by Minoo et alB 
Here m*=l/2 is the reduced mass of the electron pair. 



In Minoo et al, the physical temperature T was used and 
is valid only at high T. The use of T ee = (T 2 ±T 2 q )^ 2 
instead of T validates it down to T = as well.E£l 

Since the D + are classical particles, the D + — D + in- 
teraction is the coulomb interaction ^33 (r) = 1/r. The 
e — D + interaction is more tricky. The e-e interaction, 
eqn. ^, is based on the quantum mechanical scattering 
of two electrons. We determine the e — D + interaction 
cj> e jj(r), from the density profile n£> e {r) given by the K-S 
equation for electrons interacting with a single deuteron 
at the origin. We have discussed this in the context of 
the "neutral pseudo-atom'jJJET model (NPA-DFT) for 
solving the K-S equations .li§E3 This gives the deuteron- 
electron PDF, i.e., gDe{r) = n^, e {r)jn. This noe(r) 
includes both bound-state and continuum- state densities. 
Applying the HNC and the OZ equation to this system 
containing a single deuteron, we have, 

- (3De<pDe(r) = log[g De (r)] - h De (r) + c De {r) (7) 
hn e (r) = c De (r) + n J dr'c De (f- r')h ee (r') (8) 

The deuteron-deuteron DCF does not appear as there is 
only a single D + . Hence, knowing the gDe{r) from the so- 
lution of the Kohn-Sham equation for the single deuteron 
problem, we can obtain CDe{r) in terms of hneif)- Hence 
the e-D potential can be extracted. This determines the 
product, t3De4 } De{r), while pjj e , and <f>De are not needed 
individually. However, on solving the K-S equation for 
the regime of interest, no atomic bound states are found; 
the effective ionic charge Z — rib, where is the number 
of bound electrons per nucleus, remains unity. This does 
not contradict transient bound-states found in PIMC.EU 
Hence, at least in this regime of n and T, Kohn-Sham 
NPA-DFT is not needed; We set: 
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cj> De {r)=-{l/r)[l-e 
k De = (27rmeTee) 1/2 

l//3.D e = {T DD + T ee )/2 



(9) 
(10) 
(11) 



The first equation is just the V cou with the r = value 
set to an inverse dBL for the D-e pair, koe, as in the 
e-e interaction. The dBL 1/feoe contains only the elec- 
tron contribution since the D contribution is zero for a 
classical particle. Thus only the T ee appears in kDe (the 
effective mass of the D-e pair is m e , since the deuteron 
mass Mo » m e ). The effective temperature l//?z>e 
for the D-e interaction is the mean kinetic energy of the 
pair. However, the use of the K-S procedure for the di- 
mensionless potential ^De't'Deij') becomes obligatory at 
lower densities (r s > 2) and temperatures (T < 1.0 eV) 
when the bound-state population n& becomes nonzero. 

Using the above scheme, we solve the coupled set of 
HNC equations to determine the six PDFs of the e, D + 
system. The excess free-energy F exc (r s , T) is determined 
via a coupling-constant integration, as in ref. [ [l7|| , for 
a range of values of r s and T, and converted into the 
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total free energy F(r s ,T) by adding on the ideal elec- 
tron and ion contributions Fq,Fq (see Fig. 1(a)). The 
total pressure P and the total internal energy E are ob- 
tained as usual by P — dF(r s ,T)/dV, where V is the 
volume, and E = df3F(r s ,T) / 'd/3, where 0=1/T. In the 
regime of interest, i.e, 1.8 < r s < 2.1, and 0.8 eV < T < 
15, we find that F exc {r s ,T) is approximately linear, i.e., 
F exc (r s ,T) = M(T)r s + C{T). The T-dependence of 
M(T), C(T) is quite nonlinear. Figure 1(b) shows that 
M (T) changes character near T — 4 eV. 

Our P, E results are compared with the PIMC of Mil- 
itzer and Ceperley (MC) in table I, showing good agree- 
ment for T > 2.75 eV. For lower T, our pressures are 
smaller. 

The free-energy F(r s ,T) is used to calculate the deu- 
terium Hugoniot for the initial state, (Eg, Vo, Pq), with T 
= 19.6 K and an initial density po — 0.171 g/cm 3 . The 
initial-state Eq = -15.886 eV per atom, and Po = 0. The 
CHNC result, and others are shown in Fig. 2. The CHNC 
Hugoniot, similar to PIMC, approaches SESAME at high 
temperatures. A softening of the Hugoniot around 2 
Mbar, not seen in the PIMC curve is also noted. This ap- 
pears near n and T where the interacting electron chem- 
ical potential /i e (r s ,T) passes through zero. 

We turn to the HNC equations to consider electrons 
and deuterons at two different temperatures, T e and Tjj. 
The shock is launched from an Aluminum pusher; the 
ions are initially much hotter than the electrons. The ve- 
locity measurements begin after about 3 ns in the laser 
experiments, and after a longer time in the Sandia work. 
Landau-Spitzer theory would indicate that the D + — e 
equilibriation occurs well within nanosecond timescales. 
However, the formation of coupled modes in plasmas 
with r > 1 streuigly quenches the ion-electron equilib- 
riation process.Ea AJso, experimental evidence exists for 
this point of viewH3 A compactly held screening-charge 
at each ion would act like a neutral object which, while 
having a very hot deuteron at the centre, would screen 
it from the cooler outer electrons. The effect could lower 
the electron-ion relaxation by more than an order of 
magnitude. c3£3 In lieu of a systematic energy relaxation 
analysis, here we assume that the deuteron nuclei are 
about 5 eV hotter than the nominal electron temperature 
( i.e., Td = T e + 5 eV), except at the lowest tempera- 
tures. Using T e , T D in the HNC equations as before, we 
have calculated a quasi-equilibrium F exc (r s ,T e ,Tn) and 
a shock Hugoniot (quasi-equilibrium concepts are dis- 
cussed in ref. 22). The resulting non-equilibrium Hugo- 



niot is given in Fig. 2. 

The higher To makes the deuteron-electron fluid more 
compressible. This appears counter-intuitive if one con- 
siders only the D + contribution. The quasi free energy 
F(r s ,T e ,T D ) consists of F e , F D , and F De . On setting 
Td > T e , the Fd term taken alone reduces the compress- 
ibility, but the total compressibility is increased by the 
major role of the pair-term P£> e . As seen in Fig.l (a), the 
fluid is in a regime close to the /i e = transition. Thus 
a higher Td increases To e and reduces the electron de- 



generacy even more, making it more compressible. This 
also reduces the screening and makes the electrons in- 
teract more strongly with the nuclei. When this effect 
is strong enough to offset the reduction of the compress- 
ibility from the ideal gas term of the hotter deuterons, a 
softening of the Hugoniot could result. In Fig. 2 we show 
a Hugoniot labeled NEQO where the ideal term Po was 
computed just as in the equilibrium Hugoniot, while in 
NEQ the full effect was included. Thus we see that ex- 
cept for the lowest temperatures, the contribution of the 
deuteron-electron pairs dominate. 

Our explanation of the observed laser-shock Hugoniot 
is incomplete until detailed modeling of the passage to 
equilibrium is achieved. Other factors like the planarity, 
constancy and duration of the Shockwave are also rele- 
vant. However, the present discussion strongly justifies 
the need to include equilibriation effects as well. The 
CHNC approach presented here is numerically and com- 
putationally very simple and should be a handy tool 
in such studies. In fact all the calculations presented 
here were-parried out using very modest computational 
facilities.cJ Our computer codes may be remotely ac-. 
cessed by interested researchers by visiting our website. cB 

In conclusion, we present a parameter- free calculation 
of the EOS of deuterium in the regime of densities and 
temperatures addressed by recent laser and magnetic 
shock experiments. Kohn-Sham calculations show the 
absence of atomic bound states in this regime. Hence a 
simple classical mapping of the electron quantum fluid 
was used. The calculated Hugoniot is in good agreement 
with other first-principles calculations. Calculations for 
the case with the ions hotter than the electrons by about 
5 eV are also presented. They suggest the possibility 
that the anomalous Hugoniots obtained from laser ex- 
periments could arise from the slowness of the system to 
cquilibriate within the experimental time scales. 
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TABLE I. The total pressure P (Mbar) and total energy 
E (eV), calculated using the Classical-map HNC (CHNC) ap- 
proach, and the path-integral Monte Carlo (PIMC) approach 
of Militzer and Ceperley (MC) at r s =2.0, i.e, at a deuterium 
density of 0.6691 g/cm 3 . 



T{K) 




P(CHNC) 


P(MC) 


^(CHNC) 


E(MC) 


500 000 


-5.35310 


26.278 


25.980 


113.30 


113.20 


250 000 


-2.14960 


12.244 


12.120 


47.57 


45.70 


125 000 


-0.99712 


5.374 


5.290 


13.60 


11.50 


65 000 


-0.64405 


2.143 


2.280 


-3.21 


-3.80 


31 250 


-0.57058 


0.754 


1.110 


-10.74 


-9.90 


15 625 


-0.57119 


0.213 


0.540 


-13.97 


-12.90 


10 000 


-0.57890 


0.058 


0.470 


-14.91 


-13.60 
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FIG. 1. (a) The temperature dependence of the excess 
free energy F exc (r s ,T), F e ((r s ,T), and the interacting chem- 
ical potential fj, e (r a ,T), at r s = 1.86. Note the passage of fi 
through zero, (b) The excess free energy in the range of 
the the shock experiments can be fitted to the linear form 
F exc (r s ,T) = M(T)r 3 + C(T); the slope M(T) and the inter- 
cept C(T) are shown as a function of T. Nore the change of 
character in M(T) when jj, e passes through zero. 



FIG. 2. Comparison of the CHNC Hugoniot with experi- 
ment and other theories. Two non-equilibrium Huoniots are 
also shown (see text). Experiments 1, 2 and 3 refer to Da 
Silva et al., Collins et al, and Knudson et al, respectively. 
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